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ABSTRACT 

The two-dimensional spatially elliptic Navier-Stokes equations have been used to investigate 
the radiative interactions in chemically reacting compressible flows of premixed hydrogen and air 
in an expanding nozzle. The radiative heat transfer term in the energy equation is simulated using 
the Monte Carlo method (MCM). The nongray model employed is based on the statistical narrow 
band model with an exponential-tailed inverse intensity distribution. The spectral correlation has 
been considered in the Monte Carlo formulations. Results obtained demonstrate that the effect 
of radiation on the flowfield is minimal but its effect on the wall heat transfer is significant. 
Extensive parametric studies are conducted to investigate the effects of equivalence ratio, wall 
temperature, inlet flow temperature, and the nozzle size on the radiative and conductive wall 
fluxes. 

NOMENCLATURE 


Latin Symbols 

A reaction rate constant; also area, m 2 

C concentration, kg.mole/m 3 

C p specific heat, J/(kg.K) 

D diffusion coefficient, m 2 /s 

E total internal energy, J/kg; also activation energy, J/kg 

f mass fraction 

g Gibbs energy, J/(kg.K) 

Graduate Research Assistant. Student Member ASME. 

Eminent Professor. Fellow ASME. 
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h 

h R 

k 


kb 

ket 

k f 

L 


mu 

M 

N 

N s 

N r 

P 

Qcw 

-V.q r 

Qrw 

Q 

R 

R u 

s, s',s" 
t 

T 


u 

u 

V 


w 

X 

X 


static enthalpy, J/kg 
base enthalpy, J/kg 

spectral radiative intensity, kW/(m 2 .sr.cm‘ l ) 

thermal conductivity, J/(m.s.k); also line intensity to spacing ratio, cm' 
atm ' 1 

backward rate constant 
equilibrium constant 
forward rate constant 
nozzle length, m 
total number of narrow bands 
molecular weight 

temperature coefficient in reaction rate expression 

number of species 

number of reactions 

gas pressure, atm 

conductive wall flux, kW/m 2 

radiative source term, kW/m 3 

net radiative wall flux, kW/m 2 

radiative energy per unit volume, kW/m 3 

gas constant, J/(kg.K); also random number 

universal gas constant, J/(kg.K) 

position variables, m 

time, s 

absolute temperature, K 
velocity in x direction, m/s 
diffusion velocity in x direction, m/s 
velocity in y direction, m/s 
diffusion velocity in y direction, m/s 
diffusion velocity vector, m/s 

production rate of species, kg/(m 3 .s) 
x-coordinate, m. 
mole fraction 
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y 

yb 

Greek symbols 

P 

7 

6 

9 

P 

Li 

p 

a 

T 

ru 

<i> 

uj 

n 


y-coordinate, m 

half height of cross sectional area of nozzle, m 
line width to spacing ratio 

stoichiometric coefficient; also half-width of an absorption line, cm 

equivalent line spacing, cm -1 
polar angle 

dynamic viscosity, kg/(m.s) 
computational coordinates 
density, kg/m 3 
normal stress, N/m 2 
shear stress, N/m 2 
spectral transmittance 
equivalence ratio 
azimuthal angle 

wavenumber, cm -1 
solid angle 


INTRODUCTION 

There has been extensive research underway to develop hydrogen-fueled supersonic com- 
bustion ramjet (scramjet) propulsion systems for National Aero-Space Plane (NASP). A critical 
element in the design of scramjets is the detailed understanding of the complex flowfield present 
in different regions of the system over a range of operating conditions. Numerical modeling of 
the flow in various sections has shown to be a valuable tool for gaining insight into the nature 
of these flows. 

In a hypersonic propulsion system, combustion takes place at supersonic speeds to reduce 
the deceleration energy loss. The products of hydrogen-air combustion are gases such as water 
vapor and hydroxyl radical. These species are highly radiatively absorbing and emitting. Thus, 
numerical simulation must be able to correctly handle the radiation phenomena associated with 
supersonic flows. 
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The study of radiative transmission in nonisothermal and inhomogeneous nongray gaseous 
systems requires a detailed knowledge of the absorption, emission and scattering characteristics of 
the specific species under investigation. In absorbing and emitting media, an accurate nongray 
model is of vital importance in the correct formulation of the radiative flux equations. The 
line-by-line models are theoretically the most precise models to treat radiative heat transfer. 
But solutions of the line-by-line formulation require considerably large computational resources. 
Currently, it is not practical to apply the line-by-line models in most engineering problems. The 
wide band models are the simplest nongray models and are extensively used in radiative heat 
transfer analyses [1, 2]. By far the most popular wide band model is the exponential wide 
band model developed by Edwards [3]. The exponential wide band model accounts for discrete 
absorption bands and spectral correlations resulting from the high resolution structure. However, 
the spectral discretization used in this model is too wide and it does not take into account the low 
resolution correlations between intensities and transmissivities [4. 5], Also, the case of partially 
reflecting walls cannot be correctly modelled with this approach [3]. Recently, the narrow band 
models have begun to receive a lot of attention due to the strong requirement for accurate 
simulation of radiation. Some narrow band models can compare favorably to the line-by-line 
calculations [4, 6], and they are much simpler than the line-by-line models. In addition, the 
narrow band models do not have disadvantages usually encountered with the wide band models. 

Most of the existing analyses in radiative heat transfer start with the radiative transfer 
equation. Use of a narrow band model in this equation results in two types of spectral correlations 
[7]. One is the spectral correlation between the intensity and the transmittance within the medium. 
Another is the spectral correlation between the reflected component of the wall radiosity and the 
transmittance. In order to investigate the first type of spectral correlation, all the intermediate 
transmittances in each finite volume element of medium along the path the radiative energy 
travels must be calculated and stored to make a correlated calculation. In order to investigate 
the second type of spectral correlation, a series expansion of the wall radiosity is required [8, 
9]. Essentially, this series expansion is utilized along with a technique for closure of the series. 


4 



Consideration of the history of a finite number of reflections and approximating the remaining 
reflections by a closure method in the radiative transfer equation complicates the mathematical 
formulation and increases the computer time considerably. As the geometry considered becomes 
complicated, exact simulation of radiative heat transfer by most existing methods becomes very 
difficult for the cases with reflecting walls. 

The MCM is not directly based on the radiative transfer equation to simulate radiative 
heat transfer. This results in the MCM having features different from the other methods for 
narrow band analysis. When the radiative energy is transmitted in the medium, the spectral 
correlation only occurs between the transmittances of two different segments of the same path in 
the statistical relationship for determining the absorption location of a energy bundle [10]. For 
the case with reflecting walls, Monte Carlo treatment with a narrow band model is similar to that 
with a gray model, and the second type of spectral correlation occurring in other methods does 
not exist. If the effect of scattering is included, a new type of spectral correlation occurs in the 
scattering term of the radiative transfer equation. Treatment of this spectral correlation will be far 
more complicated than the second type of spectral correlation mentioned earlier. In such cases, 
it has been indicated that only MCM can account for scattering in a correlated manner [11]. 

The objective of present study is to apply the Monte Carlo formulations with a narrow band 
model to investigate the effects of radiation on multi-dimensional chemically reacting supersonic 
flows. Only a limited number of studies are available to investigate the interaction of radiation 
heat transfer in chemically reacting viscous and supersonic flows of molecular species. Mani 
and Tiwari [12] are the first to take into account the effects of radiation in chemically reacting 
supersonic flows. This work has been extended to include relatively more advanced chemistry 
models by Tiwari et al. [13]. In both of these studies, a tangent slab approximation was 
employed with a gray gas model. This approximation treats the gas layer as a one-dimensional 
slab in evaluation of the radiative flux. Obviously, it is impossible to obtain reliable quantitative 
predictions of radiation from this treatment. In this study, one of the most accurate nongray 
models available is employed and multi-dimensional radiative heat transfer is simulated exactly 


5 



using the MCM; the results of radiative flux are then incorporated in the Navier-Stokes equations. 
This procedure provides a more accurate prediction of the radiative effects than the previous 
studies. 

GENERAL FORMULATION 


Governing Equations 


The physical model considered is a supersonic flow of premixed hydrogen and air in an 
expanding nozzle (Fig.l). The nozzle wall is defined, as noted, by a shifted sinusoidal curve. 
The inlet temperatures of hydrogen and air are considerably high so that the chemical reaction 
takes place in the entire flowfield. The products of hydrogen-air combustion include water vapor 
and hydroxyl radical. These species are highly absorbing and emitting. To simulate the flowfield 
accurately, all important phenomena such as chemistry, radiation and turbulence should be taken 
into account. In this study, the two-dimensional nozzle flow considered is described by the Navier 
Stokes and species continuity equations which take the form in the physical coordinates as 


dU_ 8F_ dG_ H 
dt dx + dy 


where vectors U, F, G and H are represented by 
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The other terms appearing in vectors F, G, and H are defined as 
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where A = — | fi y ft = fn + fi t and k = kj 4- k t . In this study the molecular viscosity fn and 
molecular thermal conductivity kj are evaluated form the Sutherland’s law [14]; the turbulent 
viscosity fit is evaluated from the Baldwin-Lomax model and the turbulent thermal conductivity 
k t is calculated from the turbulent Prandtl number. 

In Eq. (1), only (N s — 1) species equations need to be considered since the mass fraction of 
the species is prescribed by satisfying the constraint equation 

Ns 

^ /.' = 1 (14) 

1=1 

The diffusion velocity of the ith species is obtained by solving the Stefan-Max well equation 
[15], neglecting the body force and thermal diffusion effects, as 



where D{j — D\ ■ + D\-. The molecular diffusion coefficient D\- is obtained from the kinetic 
theory [15] and the turbulent diffusion coefficient D\ ■ is evaluated from the turbulent Schmidt 
number. Equation (15) has to be applied only to (N s — 1) species. The diffusion velocity for 

Ns „ 

the remaining species is prescribed by satisfying the constraint equation /jV* = 0, which 

»=i 


ensures the consistency. 

In the energy equation, it is noted that the radiative source term — V.<? r has been moved to 
the right hand side and its treatment will be different from other terms. The simulation of this 
source term will be explained in detail later. 


Thermodynamic and Chemistry Models 


The specific heat of individual species C Pi is defined by a fourth-order polynomial in 




temperature. 




The values of the coefficients appearing in Eq. (16) are found in [16]. Knowing the specific 
heat of each species, the enthalpy of each species can be found from Eq. (12) and the total 
internal energy is computed from Eq. (11). 

Chemical reaction rate expressions are usually determined by summing the contributions 
from each relevant reaction path to obtain the total rate of change of each species. Each path is 
governed by a law of mass action expression in which the rate constants can be determined from 
a temperature dependent Arrehenius expression. In vector H, the term w{ = M{C{ represents the 
net rate of production of species i in all chemical reactions and is modelled as follows: 

N a k} N, 

Y f Y ] = !,-■■ N r (17) 
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Equation (17) represents an N r step chemical reaction and Eq. (18) is the production rate for the 
ith species. The reaction constants kf } and are calculated from the following equations: 
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fj 


AjT Nj exp 



j = l,--*N r 


(19) 
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The equilibrium constant appearing in Eq. (20) is given by 
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( 24 ) 


f- = Ai(T - InT) - | '-T 2 - yT 3 - ^T 4 

-§T 5 + F, - G.T; i = 1, --- AV 


The forward rate for each reaction is determined from Eq. (19). The hydrogen-air combustion 
mechanism used in this work is from [17], but only seven species and seven reactions are 
selected for this study. The constants Aj, Nj and Ej for these reactions are listed in Table 1. The 
coefficients for the Gibb’s free energy in Eq. (24) are available in [16]. 

RADIATION TRANSFER MODEL 


The effects of radiation on the heat transfer to the nozzle flow and wall arise through the term 
—V.q r in the energy equation and the net radiative wall flux q rw . The expressions for both — V.<? r 
and q rw are very complicated integro-differential equations and they are usually treated separately 
from Eq. (1). Before applying MCM to evaluate -V.q r and q rw , temperature, concentration of 
species, and pressure in the media should be assumed. Next, the participating media and the 
surrounding walls are divided into many rectangular volume elements and surface elements (Fig. 
2). Note that the inlet and outlet surfaces of the nozzle flow are treated as pseudoblack walls 
with the same temperature as the local gas. Use of a rectangular volume element rather than 
other geometrical element can simplify the Monte Carlo simulation. However, it also introduces 
the problem that the nozzle walls do not fall on the control volume faces of the computational 
domain. In this study these curved boundaries are approximated with ladder-like lines [18] 
as shown in Fig. 2. This practice enables the modeling of complex geometries in Cartesian 
coordinates system. Errors in the approximations of these boundaries can be reduced by using 
finer computational elements. 


For an arbitrarily chosen volume element ABCD with a volume SV and an arbitrarily chosen 
surface element EF with an area 6 A (Fig. 2), the relations for — V.q r and q rw are expressed as 

_ Qv~w + Qa-sv ~ Qsv / 25 ) 

V.<7r — ct/ v * 
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Qrw — 


Qv-sa + Qa-sa - Qsa 
8A 


(26) 


Here, Qvsv and Qvsa are the total radiant energy form the entire gas that are absorbed by 
the volume element 8V and surface element 8 A, respectively; Qasv and Qa-sa are the total 
radiant energy from the bounding walls that are absorbed by 8V and 8 A, respectively; Q&v and 
Qsa are the radiant energy emitted by 8V and 8 A, respectively. 


Evaluation of the terms Qv-sv* Qa-sv* Qsv and Qv-sa in Eqs. (25) and (26) requires a 
detailed knowledge of the absorption, emission and scattering characteristics of the specific gas. 
Several models are available in the literature to represent the absorption emission characteristics 
of molecular species. An accurate nongray model employed in this study is the statistical narrow 
band model with an exponential-tailed inverse intensity distribution and the transmittance of a 
homogeneous and isothermal column of length 1 due to gas species j, averaged over [u> — (Au>/2), 
u>+(Au;/2)], is then given by [19] 


fj, ~ exp 




(27) 


where Xj represents the mole fraction of the absorbing species j; k and f3 = 271 - 7 /8 are the 
band model parameters which account for the spectral structure of the gas. The overbar symbol 
indicates that the quantity is averaged over a finite wavenumber interval Ac*/. Parameters k and 
l/<5 generated from a line-by-line calculation have been published for H 2 O and CO 2 [ 6 , 20 , 21 ]. 
The mean half-width 7 is obtained using the parameters suggested by Soufiani et al. [ 6 ]. The 
narrow band width considered is usually 25 cm” 1 . Nonisothermal and inhomogeneous media 
are treated by using the Curtis-Godson approximation [22]. 

To simulate radiative heat transfer using the MCM, the total radiant energy in a volume 
element or surface element is assumed to be composed of many energy bundles. These energy 
bundles are similar to photons in their behavior. The histories of these energy bundles are 
traced from their point of emission to their point of absorption. What happens to each of 
these bundles depends on the emissive, scattering and absorptive behavior within the medium 
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which is described by a set of statistical relationships. The detailed discussion of the MCM 
has been provided by Howell [23]. However, not all the statistical relationships given by 
Howell are applicable while using narrow band models. An important change is the necessity 
to spectrally average radiative properties within each narrow band. This results in spectrally 
correlated formulations. For the volume element ABCD shown in Fig. 2, the total emitted 
radiant energy and major statistical relationships in conjunction with a narrow band model are 
given by [24] 
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R ‘ = (33) 

where Ibu; is the Planck spectral blackbody intensity; x* is the entering location of the intensity 
from the side AD which has a length b; 9 and ^ are the polar and azimuthal angles of the intensity 
over the path s— >s’, respectively; m^ is the total number of narrow bands; R^, R x *, R$, R^ 
are random numbers which are uniformly distributed between zero and one. The statistical 
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relationships for an energy bundle emitted from a surface element are similar to those given by 
Howell [23] and they are not listed here. 

A large number of energy bundles is considered to satisfactorily represent the radiation 
emitted by a volume or surface element. The total number of energy bundles absorbed by each 
element multiplied by the energy per bundle gives the interchange of radiation among the volume 
and /or surface elements. The values of — V.<jy and q rw can then be obtained from Eqs. (25) 
and (26), respectively. Substituting —V.q T into the energy equation, Eq. (1) can be solved. 

METHOD OF SOLUTION 


Equations (1) is written in the the physical domain (x, y) and must be transformed to 
an appropriate computational domain (£, 77 ) for solution. Using an algebraic grid generation 
technique, a highly clustered grid in the physical domain (near regions where high gradients 
exist) can be obtained. In the computational domain, Eq. (1) is expressed as 


dU dF_ dG__7r 
dt dt] 


(34) 


where 


0 = C/7, F = Fy v - Gx v 
G = Gx^~ Fy 6 H = HJ 

J - x^y v - y^x v (35) 


Here J is the Jacobian of the transformation. 

The governing equation system (34) can be stiff due to the kinetic source terms contained 
in the vector H. To deal with the stiff system, the kinetic source terms are computed implicitly 
in the temporally discrete form of Eq. (34). Once the temporal discretization is performed, 
the resulting system is spatially differenced using the explicit, unsplit MacCormack predictor- 
corrector schemes. This results in a spatially and temporally discrete, simultaneous system of 
equations at each grid point. Each simultaneous system is solved, subject to initial and boundary 
conditions. At the supersonic inflow boundary, all flow quantities are specified. At the supersonic 
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outflow boundary, nonreflective boundary conditions are used. Only the upper half of the flow 
domain is computed, as the flow is assumed to be symmetric about the centerline of a two- 
dimensional nozzle. The upper boundary is treated as a solid wall. This implies a non-slip 
boundary condition. The wall temperature is given and species mass fractions and pressure are 
extrapolated from interior grid points, by assuming a non-catalytic wall as well as the boundary 
layer assumption on the pressure gradient. Symmetry boundary conditions are imposed at the 
lower boundary, that is, at the centerline. Initial conditions are obtained by specifying inlet flow 
conditions throughout the flowfield. The resulting set of equations is marched in time, until 
steady state solutions are achieved. 

The solution procedure employed in this study is summarized as following: (a) First, the 
governing equation (1) is solved without consideration of radiation in terms of the modified 
MacCormack schemes; (b) The steady solutions for temperature, concentration of species and 
pressure are then used for Monte Carlo simulation. The computed — V.^ r form the MCM is 
based on a different grid from that used for Eq. (1), Linear interpolation and extrapolation are 
employed for the transformation of —V.q r between the two grids; (c) The transformed — V.g r is 
substituted into Eq. (1), and Eq. (1) is solved again. If the differences between the new steady 
solutions and the previous steady solutions are smaller than a designated value, the computation 
ends. Otherwise, the steps (b) and (c) are repeated until solutions converge. It is noted that there 
are two levels of numerical procedures employed here which result in two different iterative 
procedures. One is the numerical procedure for solving the Eq. (1) and solutions iterated with 
time. The other is the numerical procedure for evaluating the radiative source term using the 
MCM which results in the iteration of different steady state solutions. 

RESULTS AND DISCUSSION 

Based on the theoretical and numerical analysis described earlier, a computer code has been 
developed to simulate two dimensional supersonic chemically reacting and radiating nozzle flows 
on a Cray X-MP machine. The specific goal in this study is to investigate the effects of radiation 
on the flowfield and heat flux on the nozzle wall. By referring to [25], several problems have 
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been considered. They contain four parameters: equivalence ratio of hydrogen and air, inlet 
flow temperature, wall temperature and nozzle size. Numerical solutions are obtained for a 
variety of combinations of these parameters. In each problem, flow is introduced to the nozzle 
at the same velocity of 1230 m/s and the same pressure of 1 atm. The grid size for solving the 
governing equation is 71x41. Further refinement of the grid yields little changes in the results. 
For a given radiative source distribution, the residuals of Eq. (1) are reduced by eight orders of 
magnitude in 3,000 iterations for a typical case and the steady state solutions are considered to 
have been obtained. The corresponding CPU time is about six minutes. To check the accuracy 
of computational scheme, a preliminary calculation has been carried out for chemically reacting 
nozzle flows without consideration of radiation. The results from this study show very good 
agreement with available solutions [25, 26], 

For the temperature ranges considered, the important radiating species are OH and H 2 O. But 
OH is a much less radiation participating species compared to H 2 O. In addition, the concentration 
of OH is several times less than that of H 2 O for the problems considered. So, the contribution of 
radiation from OH has been neglected in this study. For H 2 O, there are five important absorption 
bands. All these bands have been taken into account and they consist of 295 narrow bands in the 
spectral range from 150 cm -1 to 7500 cm -1 [20]. In addition, for all the problems considered, 
the nozzle wall is assumed to be gray and the wall emissivity is taken to be 0.8. 

To assure that the statistical results make sense in the Monte Carlo simulation, two require- 
ments must be met. One is the accuracy of statistical results for a given grid. The other is the 
independence of the results on a grid. In this study, the designated statistical accuracy of the 
results is defined in such a way that when the relative statistical errors of results are less than 
±5%, the probability of the results lying within these limits is greater than 95%. Independence 
of the results on a grid is considered to have been achieved when the volume element number 
in the x direction is 20 and the volume element number in the y direction varies from 10 to 
20 according to different cross-sectional height as shown in Fig. 2. For this grid, the total 
number of energy bundles had to be 5,000,000 and the required CPU time was about two hours 
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in order to meet the designated statistical accuracy in results for a typical problem. To test the 
independence of the Monte Carlo results on the grid, the same problem was investigated with a 
finer gird in which the volume element number in the x direction was increased to 30 and the 
volume element number in the y direction was doubled. To obtain the same accurate results, 
the total number of energy bundles had to increase to 15,000,000 and the corresponding CPU 
time increased to six hours. Comparing the solutions for the two different grids, it is found that 
the difference for the net radiative wall flux was never more than 2%, and the difference for the 
radiative source term was a little higher but less than 10%. If fact, the net radiative wall flux is 
the quantity we are most interested in, and its accuracy seems more important to us. 

The grid considered for Monte Carlo computations in this study is coarser than that for 
numerical solutions of the energy equation. The intermediate values of the radiative source term 
within the grid for solutions of Eq. (1) are obtained by interpolation and extrapolation. This 
should not introduce significant errors as the radiative source term is a slowly varying function 
compared to the temperature and its derivative [ 21 ]. The major CPU time consumed is in the 
Monte Carlo simulation. Fortunately, Monte Carlo subroutine only need to be called two to 
three times to obtain the converged steady state solutions. The reason for this will be explained 
later. It is believed that the computational time for Monte Carlo simulation could be reduced 
considerably if the code is vecterized and parallelized. 

It is a common knowledge that the convective heat transfer is very strong for a supersonic 
flow. So the effects of radiation may not be very important. To determine these effects 
quantitatively, a typical problem is selected in which the equivalence ratio of hydrogen and air, 
wall temperature, inlet flow temperature and the nozzle length are taken to be <£=1.0, T w =1900 
K, Ti=1900 K and L=2.0 m, respectively. Figures 3(a)-3(c) show the temperature, pressure and 
H 2 O mass fraction distributions which are essential information to analyze the effect of radiative 
heat transfer. Similar trends in results for temperature, pressure and H 2 O mass fraction are 
also observed for other cases considered. As the premixed mixture of hydrogen and air enters 
the nozzle, an exothermic chemical reaction takes place immediately, and the temperature and 
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pressure increase abruptly and reach their peaks in a region closer to the inlet location (Figs. 
3(a) and 3(b)). During this rapid change in temperature and pressure, the mass fraction of H 2 0 
also experiences a big jump from zero to a value which varies little in the rest of the flow regime 
(Fig. 3(c)). As the flow continues to move downstream, supersonic expansion plays a major 
role, and the temperature and pressure are decreased. At the same time, the chemical reaction 
proceeds but it becomes very weak. This is why there is a little change in H 2 0 mass fraction 
in the downstream region. 

Figure 4 shows the radiative source distributions at three different locations for the case 
considered in Figs. 3(a)-3(b). At the location x/L=0.1, temperature and pressure are very high 
and there is more radiant energy emitted than absorbed. Consequently, the radiative source 
distribution is higher than at locations x/L=0.5 and 0.9. The trend in results for —V.q T at the 
location x/L=0.1 is seen to be different from the results of other locations due to a decrease 
in temperature as the distance from the center line increases. The convective heat transfer 
distributions for the same locations as in Fig. 4 have been also calculated but they are not 
plotted in Fig. (4). This is because of large differences between the convective and radiative 
results; and also due to opposite signs for convective results at different locations. In most 
regions, the absolute value of the convective heat transfer is two or three orders of magnitude 
larger than the radiative source term. This situation does not change as long as the speed of the 
flow is very high. So, the effects of radiation on the flowfield are very weak for supersonic flows. 
This confirms our expectation and also answers the question that the Monte Carlo subroutine 
only needs to be called two or three times to obtain converged steady state solutions. As a matter 
of fact, a case without radiation was considered and the differences in temperature, pressure and 
H 2 0 mass fraction between the two cases were found to be less than ±1%. 

The effects of radiation on the nozzle wall flux are quite different from the flowfield. It is 
noted the radiative wall flux is dominant over the conductive wall flux. Some specific results 
obtained for radiative and convective wall fluxes are presented here. It should be noted that the 
net radiative wall flux is the weighted average of the flux quantities in x and y direction in order 
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to reduce the errors introduced by the approximation of curved wall with ladder-like lines. 

The effects of the equivalence ratio <f> on q TW and q cw are illustrated in Fig. 5. For a specific 
<j) value, q cw is seen to increase first, reach to a peak and then go down. This is compatible 
with the trend in temperature variation as seen in Fig. 3(a). Unlike q cw , q rw is seen to increase 
with distance along the nozzle. This behavior is justifiable. In this study, the inlet and outlet 
of the flow are treated as the pseudoblack walls. The oudet flow temperatures are larger than 
the inlet flow temperatures and the outlet area is also bigger than the inlet area. In addition, 
as the flow goes downstream, the cross-sectional area of the flow increases. Consequently, the 
optical length increases. These two reasons result in higher value of q rw as the distance from 
the inlet location increases. Comparing the values of q rw and q cw for each case, it is clear that 
the radiation is predominant. Even in the inlet region, q rw is more than two times higher than 
q C w • The results for three different equivalence ratios reveal different behavior for combustions 
with lean and rich mixtures. As (j> increases from 0.6 to 1.0, the flow temperature and H 2 O mass 
fraction increase by about 10% and 50% respectively, and pressure decreases by about 5%. The 
effects of these changes result in a sizable increase in the values of q rw and q cw . However, as 
$ increases from 1.0 to 1.4, the flow pressure decreases by about 5% and H 2 O mass fraction 
increases by about 15%, but the temperature shows little change. This results in only a slight 
change in the values of q TW and q cw . 

Figure 6 shows the effects of the nozzle wall temperature on q rw and q cw . The change of 
the nozzle wall temperature is found to have little influence on the combustion, and the flow 
temperature, pressure and H 2 O mass fraction remain almost the same in most regions as T w 
varies from 1500 K to 2100 K. As a result, the magnitude of the radiant energy absorbed on the 
wall is very close for the three cases with different nozzle wall temperatures. The value of q TW 
is equal to the absorbed radiant energy minus the emitted radiant energy. So q TW with higher 
wall temperature shows lower value as seen in Fig. 6. As for as q cw is concerned, except in 
the entrance region , q cw is seen to have a little change among the cases with different wall 
temperatures. This behavior is believed to be caused by the existence of turbulence. 
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The effects of the inlet flow temperature on q TW and q cw are demonstrated in Fig. 7. 
Inspection of the distribution of the q TW value among the three cases reveals a very interesting 
feature of q Tw . The values of q rw along the wall are not monotonically increased or decreased 
with Ti. The combined effects of temperature, pressure and H 2 O mass fraction in the flow on 
radiation are responsible for this behavior. It is well known that increase of temperature, pressure 
and concentration of participating medium enhances radiation. As the T, varies from 1500 K to 
1800 K and then from 1800 K to 2100 K, the flow temperature increases by about 5% while 
the pressure and H 2 O mass fraction decrease by about 10% and 15% respectively at each stage. 
An increase in temperature tries to reinforce the radiation while a decrease of pressure and H 2 O 
mass fraction tries to reduce the radiation. So there exist two driving forces which compete with 
each other to affect the radiation. As a consequence of the competition, the lowest curve for 
q TW is seen for the case with Tj= 1 800 K and the highest values are observed for the case with 
Tj= 1500 K. Unlike q rwi the values for q cw are found to increase monotonically with Tj. This 
is because the convective wall flux is only dependent on temperature. 

Finally, the effects of the nozzle size on q TW and q cw are illustrated in Fig. 8. By changing 
the nozzle length, the geometrically similar nozzles with different sizes can be obtained. As the 
nozzle length is reduced from 2.0 m to 1.0 m and then from 1.0 m to 0.5 m, the flow temperature 
and H 2 O mass fraction are decreased by about 5% while the pressure is increased by about 2% 
at each stage. The effect of an increased pressure on the radiation is overshadowed by a decrease 
in the nozzle size, temperature and H 2 O mass fraction. So, the lower values of q rw are seen in 
the figure as the nozzle length is reduced. For the smaller nozzle size, the flow temperature may 
be lower, but the derivative of temperature is actually higher. Therefore, contrary to q TWl the 
value q cw is observed to increase with a decrease in the nozzle size. The opposite trend between 
the values of q rw and q cw brings a question about the role of radiation in heat transfer on the 
nozzle wall. With a decrease of nozzle size, the differences between the values of q rw and q cw 
are reduced and the dominance of radiation is diminished. In fact, at L=0.5, the value of q cw is 
larger than the value of q rw in some parts of the nozzle wall. It is expected that the radiation 
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will become less important and the conduction will replace the radiation as dominant mode of 
heat transfer on the nozzle wall if the nozzle size continues to reduce. 

CONCLUSIONS 

The radiative interactions have been investigated for chemically reacting supersonic flows 
of premixed hydrogen and air in an expanding nozzle. The MCM has been found to be very 
convenient and reliable tool to analyze radiative heat transfer in multi-dimensional nongray 
systems. For the chemically reacting supersonic flows, the effects of radiation on the flowfield 
can be neglected but the radiative effects on the heat transfer on the nozzle wall are significant. 
The extensive parametric studies on the radiative and conductive wall fluxes have demonstrated 
that the magnitude of the radiative and conductive wall fluxes are very sensitive to the equivalence 
ratio when the equivalence ratio is less than 1.0 but they may not be so when the equivalence ratio 
is higher than 1.0. The change in the wall temperature has little effect on the combustion. Thus, 
the radiative wall flux is decreased with an increase of wall temperature. But the conductive 
wall flux seems insensitive to the change of wall temperature. The radiative wall flux does not 
change monotonically with inlet flow temperature. Lower inlet flow temperature may yield higher 
radiative wall flux. The conductive wall flux, however, increases with an increase in the inlet 
flow temperature. The radiative wall flux decreases but the conductive wall flux increases with 
a reduction of nozzle size. For large size of nozzles, the radiative wall flux is dominant over the 
conductive wall flux. However, the situation may be reversed when the nozzle size is reduced. 
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Table 1. Hydro gen- Air Combustion Mechanism (7 species, 7 reactions) 


No. 

Reaction 

A 

N 

E 

1 

H 2 + 0 2 -> OH + OH 

1.70E+13 

0.0 

24233 

2 

H + 0 2 — > OH + 0 

1.42E+14 

0.0 

8250 

3 

OH + H 2 — ► H 2 0 + H 

3.16E+07 

1.8 

1525 

4 

0 + H 2 — ► OH + H 

2.07E+14 

0.0 

6920 

5 | 

OH + OH — * H 2 0 + 0 

5.50E+13 

0.0 

3523 | 

6 

H + OH + M — > H 2 0 + M 

2.21E+22 

-2.0 

0 | 

7 

H + H + M — * H 2 + M 

6.53E+17 

-1.0 

0 


25 










































graphical 1.13: Sal Jin 29 I&S&04 1994 



Fig. 2 Grid mesh for radiation simulation 
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ABSTRACT 

The two-dimensional spatially elliptic Navier-Stokes 
equations have been used to investigate the radiative in- 
teractions in chemically reacting compressible flows of 
premixed hydrogen and air in an expanding nozzle. The 
radiative heat transfer term in the energy equation is 
simulated using the Monte Carlo method (MCM). The 
nongray model employed is based on the statistical nar- 
row band model with an exponential-tailed inverse in- 
tensity distribution. The spectral correlation has been 
considered in the Monte Carlo formulations. Results ob- 
tained demonstrate that the radiative effects on the flow- 
fleld are minimal but radiative effects on the wall heat 
transfer are significant. Extensive parametric studies are 
conducted to investigate the effects of equivalence ratio, 
wall temperature, inlet flow temperature, and the nozzle 
size on the radiative and conductive wall fluxes. 

NOMENCLATURE 

Latin Symbols 

A reaction rate constant; also area, m 2 

C concentration, kg.mole/m 3 

C p specific heat, J/(kg.K) 

D diffusion coefficient, m 2 /s 

E total internal energy, J/kg; also activation 

energy, J/kg 

f mass fraction 

g Gibbs energy, J/(kg.K) 

h static enthalpy, J/kg 

h R base enthalpy, J/kg 

Iu, spectral radiative intensity, 

kW/(m 2 .sr.cm* 1 ) 

k thermal conductivity, J/(m.s.k); also line 

intensity to spacing ratio, cm' 1 , atm* 1 

kb backward rate constant 

k^ equilibrium constant 

k f forward rate constant 
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2 Eminent Professor. Associate Fellow AIAA. 


L nozzle length, m 

rnw total number of narrow bands 

M molecular weight 

N temperature coefficient in reaction rate 

expression 

N, number of species 

N r number of reactions 

p gas pressure, atm 

qew conductive wall flux, kW/m 2 

— V.g r radiative source term, kW/m 3 
q w net radiative wall flux, kW/m 2 
Q radiative energy per unit volume, kW/m 3 
R gas constant, J/(kg.K); also random number 

R u universal gas constant, J/(kg.K) 

$, s'.s" position variables, m 
t time, s 

T absolute temperature, K 

u velocity in x direction, m/s 

u diffusion velocity in x direction, m/s 

v velocity in y direction, m/s 

v diffusion velocity in y direction, m/s 

y diffusion velocity vector, m/s 

w production rate of species, kg/(m 3 .s) 

x x-coordinate, m 

X mole fraction 

y y-coordinate, m 

y b half height of cross sectional area of 

nozzle, m 

Greek symbols 

/? line width to spacing ratio 

7 stoichiometric coefficient; also half-width 

of an absorption line, cm -1 

<5 equivalent line spacing, cm -1 
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0 polar angle 

fi dynamic viscosity, kg/(m.s) 

£ t q computational coordinates 

p density, kg/m 3 

<r normal stress, N/m 2 

t shear stress, N/m 2 

spectral transmittance 
<f> equivalence ratio 

t/» azimuthal angle 

u wavenumber, cm -1 

ft solid angle 

INTRODUCTION 

There has been extensive research underway to 
develop hydrogen-fueled supersonic combustion ramjet 
(scramjet) propulsion systems for National Aero-Space 
Plane (NASP). A critical element in the design of scram- 
jets is the detailed understanding of the complex flowfield 
present in different regions of the system over a range 
of operating conditions. Numerical modeling of the flow 
in various sections has shown to be a valuable tool for 
gaining insight into the nature of these flows. 

In a hypersonic propulsion system, combustion takes 
place at supersonic speeds to reduce the deceleration en- 
ergy loss. The products of hydrogen-air combustion are 
gases such as water vapor and hydroxyl radical. These 
species are highly radiatively absorbing and emitting. 
Thus, numerical simulation must be able to correctly han- 
dle the radiation phenomena associated with supersonic 
flows. 

The study of radiative transmission in nonisothermal 
and inhomogeneous nongray gaseous systems requires a 
detailed knowledge of the absorption, emission and scat- 
tering characteristics of the specific species under inves- 
tigation. In absorbing and emitting media, an accurate 
nongray model is of vital importance in the correct for- 
mulation of the radiative flux equations. The line-by-line 
models are theoretically the most precise models to treat 
radiative heat transfer. But solutions of the line-by-line 
formulation require considerably large-computational re- 
sources. Currently, it is not practical to apply the line- 
by-line models in most engineering problems. The wide 
band models are the simplest nongray models and are 
extensively used in radiative heat transfer analyses [1, 
2]. By far the most popular wide band model is the ex- 
ponential wide band model developed by Edwards [3]. 
The exponential wide band model accounts for discrete 
absorption bands and spectral correlations resulting from 
the high resolution structure. However, the spectral dis- 
cretization used in this model is too wide and it does 


not take into account the low resolution correlations be- 
tween intensities and transmissivities [4. 5]. Also, the 
case of partially reflecting walls cannot be correctly mod- 
elled with this approach [3], Recently, the narrow band 
models have begun to receive a lot of attention due to 
the strong requirement for accurate simulation of radia- 
tion. Some narrow band models can compare favorably 
to the line-by-line calculations [4, 6], and they are much 
simpler than the line-by-line models. In addition, the 
narrow band models do not have disadvantages usually 
encountered with the wide band models. 

Most of the existing analyses in radiative heat trans- 
fer start with the radiative transfer equation. Use of a nar- 
row band model in this equation results in two types of 
spectral correlations [7]. One is the spectral correlation 
between the intensity and the transmittance within the 
medium. Another is the spectral correlation between the 
reflected component of the wall radiosity and the trans- 
mittance. In order to investigate the first type of spectral 
correlation, all the intermediate transmittances in each 
finite volume element of medium along the path the ra- 
diative energy travels must be calculated and stored to 
make a correlated calculation. In order to investigate the 
second type of spectral correlation, a series expansion of 
the wall radiosity is required [8, 9]. Essentially, this se- 
ries expansion is utilized along with a technique for clo- 
sure of the series. Consideration of the history of a finite 
number of reflections and approximating the remaining 
reflections by a closure method in the radiative transfer 
equation complicates the mathematical formulation and 
increases the computer time considerably. As the geome- 
try considered becomes complicated, exact simulation of 
radiative heat transfer by most existing methods becomes 
very difficult for the cases with reflecting walls. 

The MCM is not directly based on the radiative 
transfer equation to simulate radiative heat transfer. This 
results in the MCM having features different from the 
other methods for narrow band analysis. When the ra- 
diative energy is transmitted in the medium, the spectral 
correlation only occurs between the transmittances of two 
different segments of the same path in the statistical re- 
lationship for determining the absorption location of a 
energy bundle [10]. For the case with reflecting walls, 
Monte Carlo treatment with a narrow band model is sim- 
ilar to that with a gray model, and the second type of 
spectral correlation occurring in other methods does not 
exist If the effect of scattering is included, a new type of 
spectral correlation occurs in the scattering term of the 
radiative transfer equation. Treatment of this spectral 
correlation will be far more complicated than the second 
type of spectral correlation mentioned earlier. In such 
cases, it has been indicated that only MCM can account 
for scattering in a correlated manner [II]. 

The objective of present study is to apply the Monte 
Carlo formulations with a narrow band model to investi- 
gate the effects of radiation on multi-dimensional chem- 
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ically reacting supersonic flows. Only a limited number 
of studies are available to investigate the interaction of 
radiation heat transfer in chemically reacting viscous and 
supersonic flows of molecular species. Mani and Tiwari 
[12] are the first to take into account the effects of radia- 
tion in chemically reacting supersonic flows. This work 
has been extended to include relatively more advanced 
chemistry models by Tiwari et al. [13]. In both of these 
studies, a tangent slab approximation was employed with 
a gray gas model. This approximation treats the gas layer 
as a one-dimensional slab in evaluation of the radiative 
flux. Obviously, it is impossible to obtain reliable quanti- 
tative predictions of radiation from this treatment. In this 
study, one of the most accurate nongray models available 
is employed and multi-dimensionai radiative heat trans- 
fer is simulated using the MCM; the results of radiative 
flux are then incorporated in the Navier-Stokes equations. 
This procedure provides a more accurate prediction of the 
radiative effects than the previous studies. 

GENERAL FORMULATION 
Governing Equations 
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The other terms appearing in vectors F, G, and H are 
defined as 



The physical model considered is a supersonic flow 
of premixed hydrogen and air in an expanding nozzle 
(Fig. 1). The nozzle wall is defined, as noted, by a shifted 
sinusoidal curve. The inlet temperatures of hydrogen 
and air are considerably high so that the chemical re- 
action takes place in the entire flowfield. The products 
of hydrogen-air combustion include water vapor and hy- 
droxyl radical. These species are highly absorbing and 
emitting. To simulate the flowfield accurately, all impor- 
tant phenomena such as chemistry, radiation and turbu- 
lence should be taken into account. In this study, the 
two-dimensional nozzle flow considered is described by 
the Navier Stokes and species continuity equations which 
take the form in the physical coordinates as 
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where vectors U, F, G and H are represented by 
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where A = — §p, p = m + pt and k — kj - f k t . In 
this study the molecular viscosity pi and molecular ther- 
mal .conductivity -are. evaluated form the Sutherland's 
law [14]; the turbulent viscosity p t is evaluated from the 
Baldwin-Lomax model and the turbulent thermal conduc- 
tivity k t is calculated from the turbulent Prandtl number. 

In Eqs. (1), only (N, — 1) species equations need to 
be considered since the mass fraction of the species is 
prescribed by satisfying the constraint equation 
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The diffusion velocity of the ith species is obtained by 
solving the Stefan-Max well equation [15], neglecting the 
body force and thermal diffusion effects, as 



(15) 


where Dij — D\j + D\j . The molecular diffusion coef- 
ficient Djj is obtained from the kinetic theory [15] and 
the turbulent diffusion coefficient Djj is evaluated from 
the turbulent Schmidt number. Equation (15) has to be 
applied only to (N g — 1) species. The diffusion veloc- 
ity for the remaining species is prescribed by satisfying 
n. _ 

the constraint equation £ /« K - = 0, which ensures the 
1=1 

consistency. 

In the energy equation, it is noted that the radiative 
source term — V.g r has been moved to the right hand 
side and its treatment will be different from other terms. 
The simulation of this source term will be explained in 
detail later. 


Thermodynamic and Chemistry Models 

The specific heat of individual species C Pi is defined 
by a fourth-order polynomial in temperature. 


% = Ai + BiT + CiT 2 + DiT 3 + E t T* (16) 
R 


The values of the coefficients appearing in Eq. (16) are 
found in [16]. Knowing the specific heat of each species, 
the enthalpy of each species can be found from Eq. (12) 
and the total internal energy is computed from Eq. (11). 

Chemical reaction rate expressions are usually deter- 
mined by summing the contributions from each relevant 
reaction path to obtain the total rate of change of each 
species. Each path is governed by a law of mass action 
expression in which the rate constants can be determined 
from a temperature dependent Arrehenius expression. In 
vector H, the term iu< = A/,C, represents the net rate of 
production of species i in all chemical reactions and is 
modelled as follows: 


N. 


*/, ^ 


E^/Q j = 1 ,-Nr (17) 

1 = 1 i=l 


W{ — A f,Cj 

= Mi £ ( t " - Yij) 

i = 1 

,* "■ 

k u nci'- n 

m=l m = l 


(18) 


Equation (17) represents an N r step chemical reaction 
and Eq. (18) is the production rate for the ith species. 


The reaction constants kj. and h. are calculated from 
the following equations: 

k u = A i TN ‘ tz p(~i^] • 7 = 1 ,"^ ( 19 ) 

it,. = k ti /k. v ; j= 1, — K (20) 

The equilibrium constant appearing in Eq. (20) is given 
by 



( 21 ) 


where 


N. N. 

A"i=E70-E7!j; i = l, (22) 

i=l i=l 


N, N, 

AG ^ = E7y«-E7 7 = 1, -(v. (23) 
«=1 «=1 


#- = A(T-/nT)-^r 2 - 


Bi 


Rt 


Ei„ 5 


~w +Fi ~ GiT; 


Q rp 3 _ Dj jA 
6 12 

i=l, — N r (24) 


The forward rate for each reaction is determined 
fromEq. (19). The hydrogen-air combustion mechanism 
used in this work is from [17], but only seven species and 
seven reactions are selected for this study. The constants 
Aj, Nj and Ej for these reactions are listed in Table 1. 
The species Gibb’s free energy expression Eq. (24) is 
obtained from the integrations of the specific heat C Pi 
and the coefficients in Eq. (24) are available in the same 
way as in Eq. (16). 

RADIATION TRANSFER MODEL 

The radiative effects on the nozzle flowfield arise 
through the term— V.g r in the energy equation and the 
radiative effects on the heat transfer on the nozzle walls 
arise through the termg rw . The expressions for both 
— V.g r and q rw are very complicated integro-differential 
equations and they are usually treated separately from the 
governing equations. Before applying MCM to evaluate 
— V.g r and q r tu* temperature, concentration of species, 
and pressure in the media should be assumed. Next, the 
participating media and the surrounding walls are divided 
into many quadrilateral volume elements and surface 
elements (Fig. 2(a)). Note that the inlet and outlet 
surfaces of the nozzle flow are treated as pseudoblack 
walls with the same temperature as the local gases. 

For an arbitrarily chosen volume element with a 
volume 6V and an arbitrarily chosen surface element 
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with an area SA in Fig. 2(a), the relations for — V.g r 
and q rw are expressed as 


-V.g r 


Qv-tv + Qa-sv — Qsv 
SV 


(25) 


_ Qv- 6A +Qa -SA ~ Qsa 
q rw = — (26) 

Here, Qv-tv and Qv-sa are the total radiant energy 
form the entire gas that are absorbed by the volume ele- 
ment SV and surface element SA, respectively; Qa-sv 
and Qa-sa are the total radiant energy from the bound- 
ing walls that are absorbed by SV and SA, respectively; 
Qsv and Qsa are the radiant energy emitted by SV and 
SA, respectively. 

Evaluation of the terms Qv-tv, Qa-sv, Qsv and 
Qv-sa in Eqs. (25) and (26) requires a detailed knowl- 
edge of the absorption, emission and scattering character- 
istics of the specific gas. Several models are available in 
the literature to represent the absorption emission char- 
acteristics of molecular species. An accurate nongray 
model employed in this study is the statistical narrow 
band model with an exponential-tailed inverse intensity 
distribution. The transmittance predicted by this model 
in a homogeneous and isothermal column of length 1 due 
to gas species j, averaged over [w — (Au>/2), u/+(Aw/2)], 
is expressed as [19J 



/ | 2 xXrfk 'j 



(27) 


where Xj represents the mole fraction of the absorbing 
species j; k and 0 = 2iry/6 are the band model param- 
eters which account for the spectral structure of the gas. 
The overbar symbol indicates that the quantity is aver- 
aged over a finite wavenumber interval Au. Parameters 
jb and 1/6 generated from a line-by-line calculation have 
been published for H2O, CO2, CO, OH, NO, and other 
species [6, 20, 21]. The mean half-width 7 is obtained 
using the parameters suggested by Soufiani et al. [6]. 
The narrow band width considered is usually 25 cm -1 . 
Nonisothermal and inhomogeneous media are treated by 
using the Curtis-Godson approximation [22]. 

To simulate radiative heat transfer using the MCM, 
the total radiant energy in a volume element or surface 
element is assumed to be composed -of-raany energy bun- 
dles. These energy bundles are similar to photons in 
their behavior. The histories of these energy bundles 
are traced from their point of emission to their point of 
absorption. What happens to each of these bundles de- 
pends on the emissive, scattering and absorptive behavior 
within the medium which is described by a set of statis- 
tical relationships. The detailed discussion of the MCM 
has been provided by Howell [23]. However, not all the 
statistical relationships given by Howell are applicable 
while using narrow band models. An important change 


is the necessity to spectrally average radiative proper- 
ties within each narrow band. This results in spectrally 
correlated formulations. For a volume element, the total 
emitted radiant energy and major statistical relationships 
in conjunction with a narrow band model are given by 
[24] 

fHy 

Qiv =4jt£ {TQfI^Au, k )dV (28) 

i=l 


Rut — 


4tt Au k )dV 

fc=i 

Qdv 


(w n_l < (J < u n ) 

(29) 


R* = 


1 — cos 9 
2 


(30) 


ib 


(31) 

/ 

'dTZ\ 

(32) 

in T w (Lm ) V 

v dl J 


Here, 7 is the mean absorption coefficient over the kth 
narrow band and is obtained as [25] 


^ In r w k(L m ) 


(33) 


In the above equations, Lm is the mean beam length of the 
volume element; I^t is the Planck spectral blackbody 
intensity for the kth narrow band; 9 and are the 
polar and azimuthal angles of emission direction of an 
energy bundle, respectively; mo; is the total number 
of narrow bands; R^, R$, Rq,, Ri are random numbers 
which are uniformly distributed between zero and one. 
The statistical relationships for an energy bundle emitted 
from a surface element are similar to those given by 
Howell [23] and they are not listed here. 

A large number of energy bundles is considered to 
satisfactorily represent the radiation emitted by a volume 
or surface element. The total number of energy bundles 
absorbed by each element multiplied by the energy per 
bundle gives the interchange of radiation among the vol- 
ume and /or surface elements. The values of — V.g r and 
q rw can then be obtained from Eqs. (25) and (26), re- 
spectively. Substituting —V.q r into the energy equation, 
Eqs. (1) can be solved. 


METHOD OF SOLUTION 
Equations (1) is written in the the physical domain 
(x, y) and must be transformed to an appropriate com- 
putational domain (£, tj) for solution. Using an algebraic 
grid generation technique, a highly clustered grid in the 
physical domain (near regions where high gradients exist) 
can be obtained as seen in Fig. 2(b). In the computa- 
tional domain, Eqs. (1) are expressed as 



dF 8G - 
di + dr, ~ 


( 34 ) 
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where 

0 = UJ, F = Fy v - Gx,j 
G = Gx r -Fyt, H = HJ 

J = x i y n -y i x n (35) 

Here J is the Jacobian of the transformation. 

The governing equation system (34) can be stiff due 
to the kinetic source terms contained in the vector H. To 
deal with the stiff system, the kinetic source terms are 
computed implicitly in the temporally discrete form of 
Eq. (34). Once the temporal discretization is performed, 
the resulting system is spatially differenced using the ex- 
plicit, unsplit MacCormack predictor-corrector schemes. 
This results in a spatially and temporally discrete, simul- 
taneous system of equations at each grid point Each si- 
multaneous system is solved, subject to initial and bound- 
ary conditions. At the supersonic inflow boundary, all 
flow quantities are specified. At the supersonic outflow 
boundary, nonreflective boundary conditions are used. 
Only the upper half of the flow domain is computed, as 
the flow is assumed to be symmetric about the center- 
line of a two-dimensional nozzle. The upper boundary is 
treated as a solid wall. This implies a non-slip boundary 
condition. The wall temperature is given and wall species 
mass fractions and pressure are extrapolated from interior 
grid points, by assuming a non-cataiytic wall as well as 
the boundary layer assumption on the pressure gradient 
Symmetry boundary conditions are imposed at the lower 
boundary, that is, at the centerline. Initial conditions are 
obtained by specifying inlet flow conditions throughout 
the flowfleld. The resulting set of equations is marched 
in time, until steady state solutions are achieved. 

With consideration of radiative heat transfer, solu- 
tion procedures employed in this study are summarized 
as following: 

(a) First, the governing equations (1) is solved with- 
out consideration of radiation in terms of the modified 
MacCormack schemes; 

(b) The steady solutions for temperature, pressure 
and species mass fractions are then used for Monte Carlo 
simulation. The computed radiative source term — V.? r 
from the MCM is based on a different grid from that 
used for Eqs. (1). Linear interpolation and extrapolation 
are employed for the trans formation of — V.g r between 
the two grids; 

(c) The transformed -V.g P is substituted into Eqs. 
(1), and Eqs. (1) is solved again. If the differences be- 
tween two consecutive steady solutions are smaller than 
a designated tolerance, the computation ends. Otherwise, 
the steps (b) and (c) are repeated until solutions converge. 

It is noted that there are two levels of numerical 
procedures employed here which result in two different 
iterative procedures. One is the numerical procedure for 
solving the Eqs. (1) and solutions iterate with time. 
The other is the numerical procedure for evaluating the 


radiative source term using the MCM which results in 
the iteration of steady state solutions. 

RESULTS AND DISCUSSION 

Based on the theoretical and numerical analysis de- 
scribed earlier, a computer code has been developed to 
simulate two-dimensional supersonic chemically react- 
ing and radiating nozzle flows on a Cray X-MP machine. 
The specific goal in this study is to investigate the effects 
of radiation on the flowfleld and heat flux on the nozzle 
wall. By referring to [26], several problems have been 
considered. They contain four parameters: equivalence 
ratio of hydrogen and air, inlet flow temperature, wall 
temperature and nozzle size. Numerical solutions are 
obtained for a variety of combinations of these parame- 
ters. In each problem, flow is introduced to the nozzle 
at the same velocity of 1230 m/s and the same pres- 
sure of 1 atm. The grid size for solving the governing 
equations is 71x41 (upper half of the nozzle). Further 
refinement of the grid yields little changes in the results. 
For a given radiative source distribution, the residuals of 
Eqs. (1) are reduced by eight orders of magnitude in 
3,000 iterations for a typical case and the steady state 
solutions are considered to have been obtained. The cor- 
responding CPU time is about six minutes. To check the 
accuracy of computational scheme, a preliminary calcu- 
lation has been carried out for chemically reacting noz- 
zle flows without consideration of radiation. The results 
from this study show very good agreement with avail- 
able solutions [26, 27]. 

For the temperature ranges considered, the impor- 
tant radiating species are OH and H 2 O. But OH is a much 
less radiation participating species compared to HjO. In 
addition, the concentration of OH is several times less 
than that of H 2 O for the problems considered. So, the 
contribution of radiation from OH has been neglected in 
this study. For H 2 O, there are five important absorption 
bands. All these bands have been taken into account and 
they consist of 295 narrow bands in the spectral range 
from 150 cm" 1 to 7500 cm" 1 [20]. In addition, for all 
the problems considered, the nozzle wall is assumed to 
be gray and the wall emissivity is taken to be 0.8. 

To assure that the statistical results make sense in 
the Monte Carlo simulation, two requirements must be 
met. One is the accuracy of statistical results for a given 
-grid. -The other is -the -independence of the results on 
a grid. In this study, the designated statistical accuracy 
of the results is defined in such a way that when the 
relative statistical errors of results are less than ±5%, 
the probability of the results lying within these limits is 
greater than 95%. Independence of the results on a grid 
is considered to have been achieved when the volume 
element number in the x direction is 20 and the volume 
element number in the y direction is 20 as shown in 
Fig. 2(a). For this grid, the total number of energy 
bundles had to be 5,000,000 and the required CPU time 
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was about one hours in order to meet the designated 
statistical accuracy in results for a typical problem. To 
test the independence of the Monte Carlo results on the 
grid, the same problem was investigated with a finer gird 
in which the volume element number in the x direction 
was increased to 30 and the volume element number 
in the y direction was doubled. To obtain the same 
accurate results, the total number of energy bundles had 
to increase to 15,000,000 and the corresponding CPU 
time increased to three hours. Comparing the solutions 
for the two different grids, it is found that the difference 
for the net radiative wall flux was never more than 2%, 
and the difference for the radiative source term was a 
little higher but less than 10%. If fact, the net radiative 
wall flux is the quantity we are most interested in, and 
its accuracy seems more important to us. 

The grid considered for Monte Carlo computations 
in this study is coarser than that for numerical solutions 
of the energy equation. The intermediate values of the 
radiative source term within the grid for solutions of Eqs. 
(1) are obtained by interpolation and extrapolation. This 
should not introduce significant errors as the radiative 
source term is a slowly varying function compared to the 
temperature and its derivative [28]. The major CPU time 
consumed is in the Monte Carlo simulation. Fortunately, 
Monte Carlo subroutine only need to be called one to 
two times to obtain the converged steady state solutions. 
The reason for this will be explained later. It is believed 
that the computational time for Monte Carlo simulation 
could be reduced considerably if the code is vecterized 
and parallelized. 

The radiative effects on the flowfield are investi- 
gated first. It is a common knowledge that the convective 
heat transfer is very strong for a supersonic flow. So the 
effects of radiation may not be very important. To de- 
termine these effects quantitatively, a typical problem is 
selected in which the equivalence ratio of hydrogen and 
air, wall temperature, inlet flow temperature and the noz- 
zle length are taken to be 0=1.0, T w =1900 K, Tj=l900 
K and L=2.0 m. The inlet species mass fractions are 
/«, = 0.0283, fo, = 0.2264, /h 7 o - 0.0, f OH = 
0.0, fo = 0,0, f H = 0.0, /*, = 0.74529. Figures 
3(a)-3(c) show the temperature, pressure and H 2 O mass 
fraction distributions. Knowing these information is es- 
sential to analyze the efTect of radiative heat transfer. As 
the premixed mixture of hydrogen and air enters the noz- 
zle, an exothermic chemical reaction takes place immedi- 
ately, and the temperature and pressure increase abruptly 
and reach their peaks in a region closer to the inlet lo- 
cation (Figs. 3(a) and 3(b)). During this rapid change in 
temperature and pressure, the mass fraction of H 2 O also 
experiences a big jump from zero to a value which varies 
little in the rest of the flow regime (Fig. 3(c)). As the 
flow continues to move downstream, supersonic expan- 
sion plays a major role, and the temperature and pressure 
are decreased. At the same time, the chemical reaction 
proceeds but it becomes very weak. This is why there is 


a little change in H 2 O mass fraction in the downstream 
region. Computation has been also conducted for other 
cases. Similar trends in results for temperature, pressure, 
and H 2 O mass fractions for all species are also observed. 

Figure 4 shows the radiative source distributions at 
three different locations for the case considered in Figs. 
3(a)-3(b), At the location x/L=0.1, temperature and pres- 
sure are very high and there is more radiant energy emit- 
ted than absorbed. Consequently, the radiative source 
distribution is higher than at locations x/L=0.5 and 0.9. 
The trend in results for — V.g r at the location x/L=0.1 
is seen to be different from the results of other locations 
due to a decrease in temperature as the distance from 
the center line increases. The convective heat transfer 
distributions for the same locations as in Fig. 4 have 
been also calculated but they are not plotted in Fig. 4. 
This is because of large differences between the convec- 
tive and radiative results; and also due to opposite signs 
for convective results at different locations. In most re- 
gions, the absolute value of the convective heat transfer 
is two or three orders of magnitude larger than the radia- 
tive source term. This situation does not change as long 
as the speed of the flow is very high. So, the effects of 
radiation on the flowfield are very weak for supersonic 
flows. This confirms our expectation and also answers 
the question that the Monte Carlo subroutine only needs 
to be called one or two times to obtain converged steady 
state solutions. As a matter of fact, a case without radi- 
ation was considered and the differences in temperature, 
pressure and H 2 O mass fraction between the two cases 
were found to be less than ±1%. 

The radiative effects on the heat transfer on the noz- 
zle walls are investigated next. Unlike the radiative ef- 
fects on the flowfield, the effects of radiation on the 
nozzle wall flux are significant comparing those from 
conduction. Following results will demonstrate relative 
importance of radiative and conductive wall fluxes and 
how they change with equivalence ratio, wall tempera- 
ture, inlet flow temperature, and nozzle size. Here, the 
conductive wall flux is defined as 

< 3 « 

where n represents normal direction of the wall. 

The effects of the equivalence ratio 0 on q rw and 
q C u, -are -Hlustrated in -Fig. 5. For a specific 0 value, 
q cw is seen to increase first, reach to a peak and then go 
down. This is compatible with the trend in temperature 
variation as seen in Fig. 3(a). Unlike q ew , q rw is 
seen to increase with distance along the nozzle. This 
behavior is justifiable. In this study, the inlet and outlet 
of the flow are treated as the pseudoblack, walls. The 
outlet flow temperatures are larger than the inlet flow 
temperatures and the outlet area is also bigger than the 
inlet area. In addition, as the flow goes downstream, the 
cross-sectional area of the flow increases. Consequently, 
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the optical length increases. These two reasons result in 
higher value of q rw as the distance from the inlet location 
increases. Comparing the values of q rw and q ew for each 
case, it is clear that the radiation is predominant Even in 
the inlet region, q rw is more than two times higher than 
q txo . The results for three different equivalence ratios 
reveal different behavior for combustions with lean and 
rich mixtures. As <j> increases from 0.6 to 1.0, the flow 
temperature and H 2 0 mass fraction increase by about 
10% and 50% respectively, and pressure decreases by 
about 5%. The effects of these changes result in a sizable 
increase in the values of q r w and qcw However, as <f> 
increases from 1.0 to 1.4, the flow pressure decreases by 
about 5% and H 2 0 mass fraction increases by about 15%, 
but the temperature shows little change. This results in 
only a slight change in the values of q rv > and q ew . 

Figure 6 shows the effects of the nozzle wall tem- 
perature on q rw and q ew . The change of the nozzle wall 
temperature is found to have little influence on the com- 
bustion, and the flow temperature, pressure and H 2 0 
mass fraction remain almost the same in most regions 
as T w varies from 1500 K to 2100 K. As a result, the 
magnitude of the radiant energy absorbed on the wall is 
very close for the three cases with different nozzle wall 
temperatures. The value of qrw is equal to the absorbed 
radiant energy minus the emitted radiant energy. So q r w 
with higher wall temperature shows lower value as seen 
in Fig. 6. As for as q ew is concerned, except in the en- 
trance region , q ew is seen to have a little change among 
the cases with different wall temperatures. This behavior 
is believed to be caused by the existence of turbulence. 

The effects of the inlet flow temperature on q rw 
and q ew are demonstrated in fug. 7. Inspection of 
the distribution of the q r w value among the three cases 
reveals a very interesting feature of g PW . The values of 
q rw along the wall are not monotonically increased or 
decreased with T*. The combined effects of temperature, 
pressure and H 2 0 mass fraction in the flow on radiation 
are responsible for this behavior. It is well known that 
increase of temperature, pressure and concentration of 
participating medium enhances radiation. As the Ti 
varies from 1500 K to 1800 K and then from 1800 K 
to 2100 K. the flow temperature increases by about 5% 
white the pressure and H 2 0 mass fraction decrease by 
about 10% and 15% respectively at each stage. An 
increase in temperature fries to reinforce the radiation 
while a decrease of pressure and H 2 0 mass fraction tries 
to reduce the radiation. So there exist two driving forces 
which compete with each other to affect the radiation. 
As a consequence of the competition, the lowest curve 
for q rw is seen for the case with Ti= 1800. K and the 
highest values are observed for the case with Tj= 1500 
K. Unlike g rw , the values for q ew are found to increase 
monotonically with T s . This is because the convective 
wall flux is only dependent on temperature. 

Finally, the effects of the nozzle size on q rw and 


q ew are illustrated in Fig. 8. By changing the nozzle 
length, the geometrically similar nozzles with different 
sizes can be obtained. As the nozzle length is reduced 
from 2.0 m to 1.0 m and then from 1.0 m to 0.5 m, the 
flow temperature and H 2 0 mass fraction are decreased 
by about 5% while the pressure is increased by about 
2% at each stage. The effect of an increased pressure 
on the radiation is overshadowed by a decrease in the 
nozzle size, temperature and H 2 0 mass fraction. So, 
the lower values of q rw arc seen in the figure as the 
nozzle length is reduced. For the smaller nozzle size, 
the flow temperature may be lower, but the derivative 
of temperature is actually higher. Therefore, contrary 
to q rw , the value q cw is observed to increase with a 
decrease in the nozzle size. The opposite trend between 
the values of q r w and q e w brings a question about the role 
of radiation in heat transfer on the nozzle wall. With 
a decrease of nozzle size, the differences between the 
values of g rw and q cw are reduced and the dominance 
of radiation is diminished. In fact, at L=0.5, the value 
of q cw is larger than the value of q rw in some parts of 
the nozzle wall. It is expected that the radiation will 
become less important and the conduction will replace 
the radiation as dominant mode of heat transfer on the 
nozzle wall if the nozzle size continues to reduce. 


CONCLUSIONS 

The radiative interactions have been investigated for 
chemically reacting supersonic flows of premixed hydro- 
gen and air in an expanding nozzle. The MCM has been 
found to be very convenient and reliable tool to analyze 
radiative heat transfer in multi-dimensional nongray sys- 
tems. For the chemically reacting supersonic flows, the 
effects of radiation on the ftowfield can be neglected but 
the radiative effects on the heat transfer on the nozzle 
wall are significant. The extensive parametric studies 
on the radiative and conductive wall fluxes have demon- 
strated that the magnitude of the radiative and conduc- 
tive wall fluxes are very sensitive to the equivalence ratio 
when the equivalence ratio is less than 1.0 but they may 
not be so when the equivalence ratio is higher than 1.0. 
The change in the wall temperature has little effect on 
the combustion. Thus, the radiative wail flux is decreased 
with an increase of wall temperature. But the conduc- 
tive wall flux seems -insensitive -to the change of wall 
temperature. The radiative wall flux does not change 
monotonically with inlet flow temperature. Lower inlet 
flow temperature may yield higher radiative wall flux. 
The conductive wall flux, however, increases with an in- 
crease in the inlet flow temperature. The radiative wall 
flux decreases but the conductive wall flux increases with 
a reduction of nozzle size. For large size of nozzles, the 
radiative wall flux is dominant over the conductive wall 
flux. However, the situation may be reversed when the 
nozzle size is reduced. 
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0.0 

6920 

5 

OH + OH — H 2 0 + O 

5.50E+13 
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-1.0 
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Fig. 2(b) Grid mesh for flowfield simulation. 
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Fig. 3(a) Temperature contours in the nozzle. 
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Fig.3(b) Pressure contours in the nozzle. 
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Fig. 3(c) H 2 0 mass fraction contours in the nozzle. 
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Fig.4 Radiativa sourca distributions at three locations. 
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